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Fourth  Order  Partial  Differential  Equations  on  General  Geometries* 


John  B.  Greer  '  Andrea  L.  Bertozzi  *  Guillermo  Sapiros 


Abstract 

We  extend  a  recently  introduced  method  for  numerically  solving  partial  differential  equations  on 
implicit  surfaces  (Bertalmfo,  Cheng,  Osher,  and  Sapiro  2001)  to  fourth  order  PDEs  including  the  Cahn- 
Hilliard  equation  and  a  lubrication  model  for  curved  surfaces.  By  representing  a  surface  in  EiV  as  the 
level  set  of  a  smooth  function,  <f>,  we  compute  the  PDE  using  only  finite  differences  on  a  standard 
Cartesian  mesh  in  H:¥ .  The  higher  order  equations  introduce  a  number  of  challenges  that  are  of  small 
concern  when  applying  this  method  to  first  and  second  order  PDEs.  Many  of  these  problems,  such 
as  time-stepping  restrictions  and  large  stencil  sizes,  are  shared  by  standard  fourth  order  equations  in 
Euclidean  domains,  but  others  are  caused  by  the  extreme  degeneracy  of  the  PDEs  that  result  from 
this  method  and  the  general  geometry.  We  approach  these  difficulties  by  applying  convexity  splitting 
methods,  ADI  schemes,  and  iterative  solvers.  We  discuss  in  detail  the  differences  between  computing 
these  fourth  order  equations  and  computing  the  first  and  second  order  PDEs  considered  in  earlier 
work.  We  explicitly  derive  schemes  for  the  linear  fourth  order  diffusion,  the  Cahn-Hilliard  equation  for 
phase  transition  in  a  binary  alloy,  and  surface  tension  driven  flows  on  complex  geometries.  Numerical 
examples  validating  our  methods  are  presented  for  these  flows  for  data  on  general  surfaces. 


1  Introduction 

Partial  differential  equations  (PDEs)  defi  ned  on  surfaces  embedded  in  HI5  arise  in  a  wide  range  of  applica¬ 
tions,  including  fliid  dynamics,  biology  (e.g.,  fbids  on  the  lungs),  materials  science  (e.g.,  ice  formation), 
electromagnetism,  image  processing  (e.g.,  images  on  manifolds  and  inverse  problems  such  as  EEG),  com¬ 
puter  graphics  (e.g.,  water  flawing  on  a  surface),  computer  aided  geometric  design  (e.g.,  special  curves 
on  surfaces),  and  pattern  formation.  The  work  in  this  paper  is  concerned  with  fourth  order  differential 
equations,  which  have  interests  in  all  the  areas  mentioned  above.  Examples  of  physical  flaws  modeled  by 
fourth  order  PDEs  include  ice  formation  [36,  37],  fluids  on  lungs  [21],  brain  warping  [33,  53],  and  design¬ 
ing  special  curves  on  surfaces  [22,  33].  In  this  paper  we  extend  the  work  in  [6]  to  these  high  order  flaws. 
We  represent  the  surface  with  arbitrary  geometry  implicitly,  as  the  level  set  of  a  smooth  function  defi  ned 
in  all  of  the  embedding  space  M3 ,  and  rewrite  the  relevant  equations  in  terms  of  Euclidean  coordinates 
and  derivatives  of  the  level  set  function  (see  Section  2).  This  method  has  been  used  for  solving  fi  rst  and 
second  order  equations  on  static,  [6,  27,  31,  32],  as  well  as  dynamic,  [2,  57],  surfaces.  In  [6],  the  authors 
introduced  the  general  framework  and  showed  how  to  solve  second  order  linear  and  nonlinear  diffusions 
and  reaction-diffusion  equations  on  implicitly  defi  ned  surfaces.  In  [27,  31],  the  authors  solved  the  Eikonal 
equation  on  surfaces  like  those  in  [6]  (while  in  the  fi  rst  paper  the  work  is  for  triangulated  surfaces,  in  the 
second  implicit  representations  are  used).  In  these  works,  static  surfaces  were  considered.  The  authors  of 

*Work  supported  by  the  National  Science  Foundation  and  the  Offi  ce  of  Naval  Research. 

3Courant  Institute  for  Mathematical  Sciences.  New  York  University,  New  York,  10012,  greer@cims.nyu.edu 

■^Department  of  Mathematics,  UCLA,  Los  Angeles,  CA  90095,  bertozzi@math.ucla.edu 

^Electrical  and  Computer  Engineering,  University  of  Minnesota,  Minneapolis,  MN  55455,  guille@ece.umn.edu 


1 


[2]  and  [57]  solve  second  order  diffusion  equations  on  interfaces  that  arc  deforming  subject  to  an  extrin¬ 
sic  fbw.  As  discussed  in  the  above  papers  in  detail,  implicit  representations  provide  a  natural  means  for 
addressing  these  fbws  on  arbitrary  surfaces. 

Solving  PDEs  and  variational  problems  with  polynomial  meshes  involves  the  non-trivial  discretization 
of  the  equations  in  general  polygonal  grids,  as  well  as  the  diffi  cult  numerical  computation  of  other  quanti¬ 
ties  like  projections  onto  the  discretized  surface  (when  computing  gradients  and  Laplacians  for  example). 
Although  the  use  of  triangulated  surfaces  is  quite  popular,  there  is  still  no  consensus  on  how  to  compute 
simple  differential  characteristics  such  as  tangents,  normals,  principal  directions,  and  curvatures.  On  the 
other  hand,  it  is  commonly  accepted  that  computing  these  objects  for  iso-surfaces  (implicit  representa¬ 
tions)  is  simpler  and  more  accurate  and  robust.  This  problem  becomes  even  more  signifi  cant  when  we  not 
only  have  to  compute  these  fi  rst  and  second  order  differential  characteristics  of  the  surface,  but  also  have  to 
use  them  to  solve  variational  problems  and  PDEs  for  data  defi  ned  on  the  surface.  Formal  analysis  of  fi  nite 
difference  schemes  on  non-Cartesian  meshes  is  a  new  area  [3,  9,  10,  46],  whereas  fi  nite  difference  schemes 
on  Cartesian  meshes  arc  better  understood.  Note  also  that  working  with  polygonal  representations  is  di¬ 
mensionality  dependent,  and  solving  these  equations  for  high  dimensional  (>  2)  surfaces  becomes  even 
more  challenging  and  signifi  cantly  less  studied.  The  work  here  developed  is  valid  for  all  dimensions  of 
interest.  The  computational  cost  of  working  with  implicit  representations  is  not  higher  than  with  meshes, 
since  all  the  work  is  performed  in  a  narrow  band  around  the  level-set(s)  of  interest. 

The  framework  of  implicit  representations  for  solving  PDEs  on  them,  as  introduced  in  [6,  32],  enables 
us  to  perform  all  the  computations  on  the  Cartesian  grid  corresponding  to  the  embedding  function.  These 
computations  arc,  nevertheless,  intrinsic  to  the  surface.  Advantages  of  using  Cartesian  grid  instead  of  a 
triangulated  mesh  include  the  availability  of  well  studied  numerical  techniques  (that  we  will  extend  in  this 
paper,  see  below)  with  accurate  error  measures,  and  the  topological  fbxibility  of  the  surface,  all  leading 
to  simple,  accurate,  robust  and  elegant  implementations.  The  approach  is  general  (applicable  to  PDEs 
and  variational  problems  beyond  those  derived  in  this  paper)  and  dimensionality  independent  as  well.  We 
should  note  of  course  that  the  computational  framework  here  developed  is  only  valid  for  manifolds  which 
can  be  represented  in  implicit  form  or  as  intersection  of  implicit  forms.  As  mentioned  above,  several 
disciplines  of  sciences  have  numerous  problems  that  can  be  embedded  within  the  implicit  framework. 

In  a  number  of  applications,  surfaces  arc  already  given  in  implicit  form,  see  for  example  [43, 49],  there¬ 
fore,  the  framework  in  this  paper  is  not  only  simple  and  robust,  but  it  is  also  natural  in  those  applications. 
Moreover,  in  the  state-of-the-art  and  most  commonly  used  packages  to  obtain  3D  models  from  range  data 
and  from  segmented  volumetric  medical  images,  the  algorithms  output  an  implicit  (e.g.,  distance)  func¬ 
tion  (see  for  example  graphics.stanford.edu/projects/mich/  and  www.itk.org),  and  it  is  important  to  develop 
computational  frameworks  where  the  surface  representation  is  dictated  by  the  data  and  application,  and  not 
the  other  way  around.  On  the  other  hand,  not  all  surfaces  (manifolds)  arc  originally  represented  in  implicit 
form.  For  generic  surfaces,  we  need  to  apply  an  algorithm  that  transforms  the  given  explicit  representation 
into  an  implicit  one.  Although  this  is  still  a  very  active  area  of  research,  many  very  good  algorithms  have 
been  developed;  see  for  example  [17,  23,  29,  59].  This  translation  needs  to  be  done  only  once  for  any 
surface.  For  rendering,  the  volumetric  data  can  be  used  directly,  without  the  need  for  an  intermediate  mesh 
representation. 

Once  the  surface  is  in  implicit  form,  using  the  results  and  the  basic  “dictionary”  provided  in  [6,  32],  we 
can  translate  PDEs  and  variational  problems  based  on  intrinsic  characteristics  of  the  manifold,  into  PDEs 
and  variational  problems  that  depend  on  the  implicit  manifold  and  the  embedding  space.  This  translation 
is  done  in  a  systematic  and  generic  fashion. 

In  this  paper  we  consider  fourth  order  equations  on  static  surfaces.  Although  future  work  will  un¬ 
doubtedly  extend  these  methods  to  solve  related  equations  on  dynamically  changing  surfaces,  the  jump 
from  second  (as  in  [6,  32])  to  fourth  order  equations  on  static  surfaces  is  a  signifi  cant  computational  chal- 
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lenge  unto  itself.  The  computation  of  fourth  order  diffusions  in  fht  Euclidean  space  is  far  less  understood 
than  computation  of  second  order  diffusions.  The  literature  on  numerics  for  fourth  order  PDEs  is  an  active 
area  of  research  [5,  11,  19,  20,  51,  56,  60].  Some  work  on  solving  equations  with  fourth  order  diffusions 
on  surfaces  has  been  done  for  particular  examples  (see  [37]),  but  methods  for  solving  general  fourth  order 
PDEs  on  arbitrary  smooth  surfaces  remains  untouched  due  to  the  very  complicated  nature  of  these  high 
order  equations.  We  present  a  signifi  cant  initial  step  in  this  direction. 

Many  difficulties  arise  when  computing  any  fourth  order  diffusions.  The  dynamics  depend  highly 
on  the  smoothness  of  initial  data.  Boundary  conditions  are  often  difficult  to  implement;  fourth  order 
equations  require  prescribing  two  boundary  conditions  in  contrast  to  the  one  needed  for  second  order 
diffusions.  Time  stepping  is  perhaps  the  crucial  matter  for  fourth  order  diffusions.  Stability  requirements 
restrict  explicit  time  steps  for  fourth  order  diffusions  to  be  on  the  order  of  h4.  where  h  is  the  grid  size. 
Compare  that  to  second  order  diffusions,  where  time  steps  are  only  restricted  to  be  0(h2).  The  time  step 
restriction  prohibits  explicit  schemes  in  any  meaningful  simulation  of  fourth  order  equations.  On  the  other 
hand,  implicit  schemes  require  the  inversion  of  a  linear  system  of  equations  that  is  typically  very  large  for 
fourth  order  equations.  Furthermore,  the  literature  on  solving  fourth  order  diffusions  in  arbitrary  smooth 
domains  is  virtually  nonexistent.  To  date,  a  number  of  tools  have  proved  successful  for  fourth  order 
diffusions  occurring  in  areas  such  as  thin  fi  lm  fbid  fbw  and  materials  science,  including  ADI,  spectral, 
and  finite  element  methods,  as  well  as  convexity  splitting  [5,  11,  19,  20,  51,  56,  60].  We  address  each 
of  these  methods,  and  explain  how  they  come  into  play  with  our  own  problem  when  addressing  general 
geometries. 

Although  our  work  shares  many  features  with  level  set  methods  used  for  solving  surface  diffusion 
[8,  51]  and  Willmore  fbw  [15],  there  are  signifi  cant  differences.  Both  dynamics  arc  driven  by  fourth  order 
equations,  but  unlike  the  problems  we  consider  here,  in  surface  diffusion  and  Willmore  fbw  the  surface 
is  given  by  the  PDE  dynamics.  Our  equations  arc  fundamentally  different  as  the  PDE  is  in  some  sense 
given  by  the  surface.  In  the  case  of  level  set  methods  for  surface  diffusion,  the  dynamics  quickly  smooth 
the  surface  while  damping  regions  of  high  curvature.  In  our  case,  those  regions  of  high  curvature  remain 
unchanged  and  stay  an  integral  part  of  the  PDEs  that  we  solve.  Our  computational  methods  also  differ 
from  [8],  [15],  and  [51].  Unlike  the  work  in  [8],  we  use  implicit  times  stepping  schemes.  Also,  since 
we  solve  only  in  a  band  around  the  surface  to  speed-up  the  computation,  we  can  not  directly  use  the  FFT 
method  discussed  in  [51].  However  we  use  operator  splitting  schemes  much  like  those  in  [51].  The  authors 
of  [15]  use  fi  nite  element  methods,  whereas  we  use  fi  nite  difference  schemes. 

In  the  next  section  we  fi  rst  provide  the  basic  structure  of  the  problem,  explain  its  challenges,  and 
describe  the  particular'  equations  we  will  address.  The  general  organization  of  the  paper  is  detailed  after 
this. 


2  Basics,  Challenges,  and  Goals 


In  order  to  adequately  describe  the  challenge  of  computing  higher  order  equations  on  surfaces,  we  fi  rst 
summarize  the  ideas  developed  in  [2,  6,  57],  and  then  describe  where  those  methods  must  be  modified. 
Let  5  be  a  given  smooth  closed  curve  or  surface  in  R;V ,  where  N  is  assumed  to  be  two  when  S  is  a  curve, 
or  three,  when  5  is  a  surface.  Let  4» :  — >  R  be  a  smooth  function  whose  zero  level  set  is  given  by  S  with 
4>  <  0  inside  S  and  c]>  >  0  outside  S.  Thus  gives  the  outer  pointing  normal  to  .S'.  Letting  I  denote  the 
identity  matrix,  defi  ne 


V(|)<g)  V(f) 

|V<|>|2 


(1) 
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so  that  P(xo)v(xo)  gives  the  projection  of  a  vector  v  at  the  point  xo  6  M'v  onto  the  tangent  plane  to  the 
surface 

5  =  {x  E  M3  :  <|)(x)  =  4> (xo ) } • 

Suppose  u  is  a  smooth  function  defi  ned  on  S.  We  use  Vsu  to  denote  the  gradient  of  u  intrinsic  to  the 
surface  S.  This  is  simply  the  projection  onto  S  of  the  gradient  of  (the  extended)  u  in  the  embedded  space 
[50].  We  similarly  use  A $f  to  denote  the  Laplacian  of  it  intrinsic  to  .S’;  in  other  words,  Ay  is  the  Laplace- 
Beltrami  operator  on  S.  Now  suppose  u  is  a  function  defined  on  all  of  Mv .  We  use  Vw  and  Am  to  denote 
the  the  standard  Euclidean  gradient  and  Laplacian  of  u.  We  calculate  V<,  and  A<,  using  P  and  extrinsic 
derivatives  (i.e.,  derivatives  in  the  Euclidean  space  containing  S):  for  all  points  x  on  5, 


Vsu(x)  =P(x)Vm(x)  | V<|)(x) I- 


(2) 


Since  (])  and  P  are  defi  ned  in  some  domain  £2  c  i?  containing  S,  if  u  is  defi  ned  in  all  of  £2,  then  equation  (2) 
makes  sense  in  all  of  £2.  At  each  x  €  £2,  V  su(x)  defi  ned  as  above  corresponds  to  the  gradient  of  u  intrinsic 
to  the  level  set  of  (|)  through  x.  We  compute  the  Laplace-Beltrami  operator  in  the  same  manner: 


Asn(x) 

If  (f)  is  a  signed  distance  function, 
then  (2)  and  (3)  simplify  to 


I^|V-(P(x)Vm(x)  |V4>(x)|). 


m  =  if 


Vs«(x)  =P(x)Vm(x) 


and 

A su(x)  —  V  •  (P(x)Vm(x)  ) . 

It  is  thus  desirable,  though  unnecessary,  to  defi  ne  <|)  as  a  signed  distance  function. 


(3) 


2.1  Example:  The  Heat  Equation 

Consider  the  heat  equation  on  a  given  surface  S.  This  equation  arises  from  the  gradient  descent  [6]  of  the 
energy 

E{u)  =  fl\Vsufa  (4) 

which  is  given  by 


ut  =  Vs  •  V5M  =  A su. 


(5) 


We  combine  (5)  with  an  initial  condition, 

u(y,  0)  =  /(}’)  for  yon  5, 


(6) 


to  derive  the  heat  equation  on  the  surface  S.  We  assume  (|)  is  a  distance  function  as  discussed  above  and 
use  (3)  to  write  (5)  in  terms  of  Euclidean  derivatives: 


ut  —  V  •  (P  Vm)  .  (7) 

Equation  (7)  is  defi  ned  in  all  of  and  can  therefore  be  computed  on  a  Cartesian  grid.  Solving  (7)  in  all  of 
M'v  is  equivalent  to  solving  (5)  on  each  level  set  of  (|).  In  particular,  it  will  be  necessary  to  defi  ne  an  initial 
condition  no  on  the  entire  computational  domain  (a  band  around  S )  that  is  equal  to  /  on  S.  Equation  (7)  is 
discussed  and  computed  in  [2,  6,  57]. 
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2.2  Example:  Linear  Fourth  Order  Diffusion 

This  paper  centers  on  fourth  order  equations,  the  simplest  example  being  linear  fourth  order  diffusion. 


ut  =  ~AsAsu  (8) 

m(x,0)  =  «(). 

Equation  (8)  is  the  gradient  descent  of  the  energy 

E{u)  =  ^JjAsu\2. 

If  S  is  given  by  the  level  set  of  a  distance  function  (|),  then  we  can  compute  (8)  in  the  Euclidean  space 

containing  S  (and  then  use  Cartesian  numerics)  by  the  equation 

ut  =  — V  ■  (P  V  (V  ■  (P  Vw)))  (9) 

m(x,0)  =  Mo, 

where  P  is  the  projection  operator  given  by  (1). 

2.3  Basic  Observations  and  Challenges 

Equation  (9)  shares  a  number  of  challenges  with  (7).  These  are  problems  that  have  been  addressed  for  the 

second  order  problem,  but  take  on  additional  features  and  challenges  in  the  high  order  equation: 

1.  Domain.  We  first  note  a  key  element  of  equations  (7)  and  (9)  as  compared  with  equations  (5)  and 
(8).  The  new  PDEs  are  defi  ned  in  all  of  1#.  This  is  both  the  advantage  of  the  method  and  one  of  its 
main  challenges.  Since  the  new  PDEs  arc  defi  ned  in  Eulerian  coordinates,  we  can  use  a  Cartesian 
grid  and  apply  the  usual  fi  nite  difference  schemes  for  solving  these  surface  equations.  However  since 
the  PDEs  arc  defi  ned  in  all  of  1# ,  we  have  increased  the  problem  by  a  dimension.  We  minimize  the 
additional  work  by  computing  only  in  a  band  around  the  surface.  Fixing  c  >  0  and  considering  a 
signed  distance  function  c|)(x),  we  compute  only  in  the  band  |<|)(x)|  <  c.  Unfortunately,  this  requires 
computing  on  unusual  domains  with  curved  boundaries,  while  to  date,  most  simulations  of  fourth 
order  equations  have  been  done  on  rectangular  domains.  Thus  we  must  simultaneously  develop 
methods  for  solving  fourth  order  equations  on  Cartesian  grids  with  complicated  boundaries. 

We  also  need  to  choose  appropriate  initial  data  for  (7)  and  (9).  For  the  underlying  problem,  our 
initial  data  is  only  defi  ned  on  the  surface.  We  must  extend  these  values  to  a  function  defi  ned  in  the 
entire  band. 

2.  Degeneracy.  Consider  equation  (7).  At  each  point  x  £  M.N ,  (7)  defi  nes  a  diffusion  that  is  degenerate 
in  one  direction.  There  is  no  diffusion  in  the  direction  normal  to  A,  so  the  equation  is  extremely 
anisotropic.  Similar-  anisotropic  second  order  diffusions  have  been  thoroughly  studied;  consider  for 
example  anisotropic  diffusions  used  for  image  processing.  One  might  argue  that  the  fourth  order 
equation  (9)  is  even  more  degenerate,  since  it  is  a  higher  order  diffusion  with  absolutely  no  diffusion 
in  one  direction.  To  our  knowledge,  no  such  equation  has  been  studied.  The  closest  example  we 
know  of  is  surface  diffusion. 

This  degeneracy  plays  a  central  role  in  our  choice  of  computational  methods.  It  leads  us  to  consider 
convexity  splitting  methods  similar-  to  those  used  in  [57]  to  deal  with  the  degeneracy  of  the  second 
order  problem.  In  Appendix  B  we  look  at  an  ADI  method  that  has  been  previously  suggested  for 
fourth  order  diffusions  [56].  We  examine  this  scheme  and  show  why  the  degeneracy  of  the  fourth 
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order  problem  prohibits  a  direct  application  of  this  method.  We  also  show  how  the  second  order 
degeneracy  is  better  behaved,  as  it  is  amenable  to  the  same  approach.  Although  ADI  can  not  be 
applied  directly  to  equation  (9),  we  do  use  it  to  invert  the  linear  biharmonic  operator  in  convexity 
splitting  schemes  that  we  discuss  in  Section  4. 

3.  Boundary  Conditions.  Suppose  we  arc  solving  on  the  domain  Sc  =  {x|(|)(x)  <  c} .  No  boundary 
conditions  arc  required  to  solve  the  above  PDEs  in  Sc,  since  no  information  is  shared  between  level 
sets.  The  values  of  u  on  the  boundaries  (j>  =  ±c  thus  have  no  effect  on  u  in  the  interior  of  the 
band.  However,  since  we  discretize  on  a  Cartesian  grid,  any  numerical  solution  will  depend  on  the 
boundary,  though  this  dependence  should  decrease  as  It.  At  — >■  0.  We  also  note  that  solving  fourth 
order  PDEs  on  bounded  domains  requires  prescribing  two  boundary  conditions;  for  example  one 
might  fi  xu  and  A u,  or  Vw  •  v  and  V Ait  ■  v  on  the  boundary  (with  normal  v). 

Although  we  build  upon  previous  work  on  equations  like  (7),  we  do  not  simply  re-apply  those  meth¬ 
ods.  Equation  (9)  inherits  all  of  the  complications  of  other  fourth  order  equations  and  has  a  few  that  are 
particular  to  its  nature  as  a  higher  order  surface  PDE.  These  complications  have  not  been  addressed  in 
the  work  on  second  order  surface  PDEs,  and  in  fact  many  of  them  are  central  to  current  research  on  the 
computation  of  general  higher  order  equations: 

1.  Dependence  on  high  derivatives.  This  obvious  difference  between  equations  (7)  and  (9)  presents 
many  diffi  culties  that  were  not  apparent  in  the  works  discussed  above.  Equation  (9)  depends  on 
fourth  derivatives  of  not  only  n,  but  of  4>  as  well.  This  high  order  dependence  places  severe  re¬ 
strictions  on  u  and  (f).  Both  must  be  smooth  and  accurate  to  a  high  enough  degree  that  they  possess 
smooth  numerical  derivatives.  As  we  shall  see,  level  set  functions  cf)  that  are  smooth  enough  to  use  in 
computations  of  second  order  PDEs  might  not  be  smooth  enough  for  their  fourth  order  counterparts. 
The  initial  condition  uq  must  also  be  suffi  ciently  smooth.  This  plays  a  particularly  important  role 
when  extending  the  initial  condition  on  S  to  a  band  containing  S.  See  section  3.4  for  more  details. 

2.  Stencil  Size.  In  one  dimension,  the  standard  centered  difference  discretization  of  the  biharmonic 
operator  has  a  stencil  width  of  5  grid  points.  In  three  dimensions,  the  stencil  size  is  33  grid  points. 
Compare  that  to  the  7  grid  points  in  the  standard  discretization  of  the  three  dimensional  Laplacian. 
The  stencil  size  becomes  especially  problematic  when  solving  PDEs  on  a  surface.  In  all  likelihood, 
all  but  one  or  two  of  these  grid  points  will  lie  on  neighboring  contours  of  (j)  ((f)  ^2  0),  and  not  on  the 
surface  of  interest  (()  =  0).  Thus  we  must  extend  the  initial  data  off  of  S  carefully.  We  pick  uq  to 
be  constant  in  directions  perpendicular  to  S  for  this  reason.  Choosing  u  to  be  non-constant  in  this 
direction  could  signifi  cantly  increase  the  discretization  error. 

Stencil  size  also  plays  a  central  role  in  computation  speed.  Discretizing  the  Laplacian  requires  far 
fewer  operations  than  discretizing  the  biharmonic.  There  is  an  even  greater  difference  in  the  number 
of  operations  required  between  (7)  and  (9).  Computing  (7)  requires  at  least  one  matrix  multiplication 
by  P  at  each  time  step,  whereas  equation  (9)  requires  at  least  two  such  multiplications.  We  must  take 
this  need  for  heavy  computation  into  consideration  when  considering  the  domain  of  computation, 
much  more  so  than  in  the  case  of  (7). 

3.  Time  stepping.  Explicit  methods  for  fourth  order  diffusions  have  a  time  step  restriction  requiring 
the  time  step  k  to  be  on  the  order  of  h4 .  where  h  denotes  the  grid  cell  width.  This  restriction  is  far 
worse  than  the  stability  requirement  for  second  order  diffusions.  Due  to  this  drawback  of  explicit 
schemes,  we  turn  to  implicit  schemes  that  require  the  inversion  of  a  linear  system  at  each  time 
step.  Since  fourth  order  equations  depend  on  high  derivatives,  we  desire  a  highly  resolved  grid  for 
accuracy,  however  refi  ned  grids  are  detrimental  to  time  stepping.  Explicit  schemes  would  be  too 
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slow,  while  the  necessary  inversion  for  implicit  schemes  might  be  too  computationally  expensive. 
Time  stepping  remains  a  major  challenge  for  solving  fourth  order  equations. 

2.4  Overview  of  computational  examples 

We  consider  three  fundamental  fourth  order  problems  in  this  paper:  linear  fourth  order  diffusion, 

ut  =  -AsAsu, 

the  Cahn-Hilliard  equation  for  phase  transitions  in  a  binary  alloy, 

ut  —  A s  (— £2As  u  +  u3  —  u)  , 

and  a  PDE  derived  in  [47]  for  the  motion  of  thin  ti  1ms  driven  on  a  surface  by  gravity  and  surface  tension, 
§  =  n^V5K-iu4(Kl-K)-VK  (10) 

BoVs  ■  u3 gs  -  w4  ^Kl  +  •  gs  +  gnu3Vsu  ■ 

In  equation  (10),  k\  and  ki  are  the  principle  curvatures  of  the  substrate,  K  =  k\  +  K  is  curvature  of  the 
free  surface,  C,  —  u—  \ku2  +  ^k] kou3 .  and  K  is  the  curvature  tensor  for  the  curved  substrate.  The  Bond 
number  Bo  quantili  es  the  relative  strength  of  gravity  to  surface  tension.  The  vector  gs  is  the  component  of 
gravity  tangent  to  the  surface  .S’,  and  g„  is  the  magnitude  of  the  component  of  gravity  normal  to  ,S. 

For  these  problems,  we  address  each  of  the  issues  discussed  in  Section  2.3,  describe  the  methods  we 
used  to  deal  with  those  issues,  and  mention  some  areas  of  possible  future  research.  We  compare  with 
second  order  problems  throughout  the  discussion.  In  Section  3  we  discuss  the  groundwork  for  implemen¬ 
tation,  including  data  structures,  boundary  conditions,  extension  of  the  initial  data  off  of  S',  and  intermittent 
re-extension  of  the  surface  data.  Although  these  elements  have  been  discussed  previously  in  [2,  6,  57],  they 
take  on  new  features  and  challenges  with  the  fourth  order  problem.  Sections  4  and  6  concern  the  main  chal¬ 
lenge:  time  stepping  of  the  fourth  order  highly  degenerate  diffusions.  We  discuss  both  ADI  and  convexity 
splitting  methods  that  have  been  successful  for  other  fourth  order  equations  and  consider  their  relation  to 
our  problem.  The  new  challenges  here  include  the  extreme  degeneracy  of  these  surface  equations  and  the 
unusual  computational  domains. 

The  remaining  sections  demonstrate  the  generality  of  our  methods  by  applying  them  to  each  of  the 
above  equations.  In  Section  5  we  consider  linear  fourth  order  diffusion  (9)  on  the  unit  circle  in  M2  and 
the  unit  sphere  in  M2 .  We  give  pseudocode  describing  the  basic  steps  of  the  method.  We  also  check 
convergence  in  the  case  of  fourth  order  diffusion  on  a  circle,  for  which  we  know  exact  solutions.  Section  7 
concerns  the  Cahn-Hilliard  equation.  We  give  a  numerical  scheme  and  discuss  changes  to  the  algorithm 
given  for  the  linear  problem.  These  changes  arc  minor,  since  the  highest  order  term  is  still  linear.  Finally, 
in  Section  8  we  apply  these  methods  to  a  fully  nonlinear  model  for  thin  fi  lm  Ibid  fbw  on  surfaces.  We 
consider  flaws  driven  by  gravity  and  flaws  driven  by  curvature  alone.  To  demonstrate  the  effects  caused  by 
the  curvature  terms,  we  compare  the  model  that  takes  curvature  into  consideration  with  a  simpler  model 
that  ignores  curvature  effects. 

3  General  Implementation  Issues 

We  fi  rst  describe  methods  for  the  general  implementation  of  PDEs  on  implicit  surfaces,  including  data 
structures,  extension  of  initial  data  off  of  the  surface,  boundary  conditions,  and  modifi  cations  of  the  meth¬ 
ods  in  [2,  6,  57]  that  arc  required  to  compute  higher  order  PDEs.  These  arc  the  essential  building  blocks  for 
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the  numerical  computations  in  sections  5-8.  A  number  of  these  tools  can  be  used  to  solve  many  different 
PDEs  on  a  variety  of  geometries. 

We  assume  that  we  arc  given  a  level  set  function  (f)  that  determines  S  by  S  =  {x|(|)(x)  =  0}  with  (j>(x)  >  0 
for  points  x  outside  of  S  and  (j>(x)  <  0  for  x  inside  S.  Works  such  as  [43,  45,  52]  and  the  papers  mentioned 
in  the  introduction  discuss  means  of  producing  such  an  implicit  representation  (e.g.  via  highly  accurate 
WENO  schemes).  For  fourth  order  equations,  very  high  accuracy  is  required  for  the  level  set  function 
since  we  must  take  fourth  order  derivatives  of  <|>  as  well  as  of  u.  As  noted  in  Section  2,  choosing  4»  to  be  a 
signed  distance  function  simplili  es  the  PDEs. 

3.1  Surface  Complexity 

The  examples  in  this  paper  fall  into  two  categories  of  implicit  surfaces.  The  fi  rst  are  “simple”  surfaces, 
with  small  curvature,  in  which  we  have  an  analytical  expression  for  the  level  set  function.  Such  examples 
include  ellipses,  the  unit  circle  in  M2 ,  and  the  torus  and  sphere  in  M3 .  The  second  category  are  complex 
surfaces  defi  ned  computationally  by  a  dataset,  with  regions  of  high  curvature.  Our  example  of  such  a 
surface  is  the  Stanford  Bunny  [30].  In  this  second  case,  we  found  that  the  level  set  function  from  the 
original  Stanford  Bunny  data  is  not  smooth  enough  for  fourth  order  differences.  To  address  that  problem, 
we  slightly  smoothed  the  bunny  by  applying  a  few  steps  of  the  heat  equation  to  the  level  set  function.  This 
results  in  smoothing  high  derivatives  of  (])  while  retaining  the  essential  bunny  shape.  Such  smoothing  is 
not  required  for  solving  second  order  equations,  and  might  be  avoided  for  high  order  fbws  by  using  local 
mesh  refi  nement  or  adaptive  grids. 

Surface  complexity  also  affects  the  choice  of  time  stepping  method.  The  geometry  of  the  bunny  has 
regions  of  very  high  curvature  that  are  not  completely  resolved  by  our  3D  grid  with  159  x  161  x  129  points. 
The  fi  nite  difference  discretization  of  the  projected  PDE  introduces  terms  with  third  and  fourth  derivatives 
of  (|)  that  produce  a  time  scale  in  the  calculation  on  the  order  of  A.  A  stable  method  requires  resolving 
this  time  scale,  thus  introducing  a  time  step  limitation  due  to  the  geometry  of  the  surface,  as  opposed  to 
the  underlying  scheme.  When  smaller  time-steps  are  required  by  the  surface  geometry,  we  found  ADI 
schemes  to  be  more  effi  cient  than  iterative  solvers.  We  discuss  this  further  in  sections  4  and  6. 

3.2  Computational  Domain:  A  Band  Containing  S 

To  reduce  computational  costs,  we  restrict  computation  to  a  band  around  the  surface, 

Sc  =  {|<t>(x)|  <  c},  c>  0. 

The  value  of  the  parameter  c  selected  generally  depends  on  the  size  of  the  numerical  stencil  and  the 
complexity  of  the  working  surface.  This  restriction  to  a  band  is  similar  in  spirit  to  the  local  level  set 
method  [1,  45].  For  example,  when  computing  linear  fourth  order  diffusion  on  the  unit  circle,  we  could 
compute  (9)  in  a  square  containing  the  circle  ([—2,2]  x  [—2,2]  for  example),  but  instead  it  is  more  effi  cient 

to  choose  the  annulus  _ 

1  —  c  <  \J x2  +  y2  <  c  +  1 

as  the  computational  domain  (see  Fig.  1). 

3.3  Data  Structures 

Throughout  we  rely  only  on  the  implicit  structure  of  arrays.  We  follow  the  example  of  the  local  level  set 
algorithm  described  in  [45].  We  fi  rst  defi  ne  an  array  that  stores  the  values  of  <[>  within  some  box  containing 
S.  For  our  discussion  we  assume  an  A-dimensional  box  with  M  grid-points  on  each  side,  so  (])  is  stored  in 


Figure  1 :  S  denotes  the  surface  corresponding  to  the  zero  level  set  of  <|),  and  Sc  is  a  band  around  that  surface 
used  for  the  computations. 


an  array  with  MN  elements.  We  then  make  a  list  of  all  points  that  arc  within  the  band  Sc.  For  simplicity, 
we  also  store  the  values  of  u  and  any  other  functions  necessary  for  solving  the  PDE  in  arrays  of  length 
MN .  Particular  entries  of  each  array  correspond  to  particular  grid-points  in  the  box  containing  5,  however, 
to  minimize  memory  use,  one  could  use  only  one  array  of  length  MN  and  at  grid-points  lying  in  the  band, 
store  pointers  to  the  necessary  function  values.  By  using  array  storage  we  access  each  grid-point’s  nearest 
neighbors  without  any  searching.  This  is  useful  for  PDEs  discretized  by  fi  nite  differences.  To  illustrate 
how  this  would  work  in  practice,  for  our  simple  example  of  the  annulus  in  R2,  we  defi  ne  (|)  in  all  of 
Q.  =  [—2,2]  x  [—2,2]  on  a  Cartesian  grid  with  M2  grid-points  and  cell  width  h .  =  4/(M  —  1 ) .  We  mark  the 
grid-points  x/j  satisfying  |(])(x,-j|  <  c.  For  simplicity  we  also  define  u  on  the  entire  grid,  but  we  compute 
using  only  the  values  of  u  within  the  band. 

3.4  Extension  of  Initial  Surface  Data 

Since  no  is  assumed  to  be  defi  ned  only  on  the  surface  S,  we  must  extend  no  to  the  other  level  sets  of  <]) 
within  the  band.  We  extend  no  by  requiring 

Vw0-V(|)  =  0.  (11) 

In  other  words,  we  require  no  to  be  constant  perpendicular  to  the  surface.  In  our  simple  example  of  the 
annulus,  condition  (11)  states  that  in  polar  coordinates  (r,  0),  n  is  a  function  of  0  alone. 

One  way  of  forcing  (11)  is  to  solve 


,  ,  V<|> 

u,  +  Hh  (4>)  |^|  -Vn  =  0 


(12) 


to  steady  state  within  the  band,  where  Hf,(s)  is  a  smoothed  version  of  the  Heaviside  function  depending  on 
the  gridsize  h. 


Hh{s) 


s 

vA2  +  h2 
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Since  H{ 0)  =  0,  values  on  S  remain  unchanged  while  the  information  is  propagated  both  inside  (where 
H(§)  <  0)  and  outside  (//(<]))  >  0)  the  surface.  Discussion  of  this  method,  including  the  selection  of  ///,, 
can  be  found  in  [43]  and  [45].  For  second  order  equations  on  surfaces,  equation  (12)  can  be  solved  using 
simple  fi  rst  order  up-winding.  However  this  is  inadequate  for  our  purposes,  as  such  simple  integration 
results  in  an  extended  u  that  has  very  noisy  high  derivatives.  We  compute  (12)  using  fi  fth-order  WENO 
[25 ,  26]  for  spatial  derivatives  and  integrate  in  time  with  second  order  Runge-Kutta. 

3.5  Re-extension  of  Surface  Data 

We  also  note  that  u  must  be  intermittently  re-extended  off  of  the  surface,  so  the  method  used  to  compute 
(12)  must  be  both  accurate  and  effi  dent,  further  suggesting  the  use  of  a  high  order  WENO  scheme.  The 
frequency  of  re-extension  depends  on  time  steps  used  in  the  integration  of  the  surface  PDE.  If  large  time 
steps  are  taken,  we  re-extend  each  time  step  by  solving  (12)  to  steady  state.  This  can  often  be  done  in  one 
or  two  iterations  of  the  numerical  scheme  for  (12).  When  taking  small  time  steps,  we  only  re-extend  every 
few  steps  as  needed.  This  need  to  re-extend  surface  data  is  also  discussed  in  [57],  In  the  simple  example 
of  an  annular  domain,  the  re-extension  is  a  means  of  adjusting  the  values  of  u  off  of  S  so  that  it  depends 
only  on  9.  In  general,  doing  a  re-extension  enables  us  to  minimize  effects  on  surface  values  of  u  by  values 
of  u  off  of  the  surface. 

3.6  Boundary  Conditions 

Consider,  for  example,  equation  (9)  defined  only  in  the  band  Sc  —  {x|  |(j)(x)|  <  c}.  The  boundary  of  Sc 
consists  of  the  disjoint  curves  (|)  =  ±c  (see  Figure  1).  Since  the  boundary  is  given  by  level  sets  of  ([), 
information  is  not  transferred  from  the  boundary,  and  thus  boundary  conditions  are  not  needed  to  solve  the 
PDE.  However,  when  the  PDE  (9)  is  discretized  using  fi  nite  differences  on  a  Cartesian  grid  in  Sc.  because 
the  grid  directions  do  not  run  parallel  to  the  boundary  of  the  annulus,  we  must  impose  numerical  boundary 
conditions  for  computation  which  then  affect  the  solution.  However  for  short  times  and  a  large  enough 
band  size,  these  boundary  conditions  should  have  minimal  effect  on  the  value  of  the  solution  at  the  surface 
itself.  We  have  some  choice  in  the  selection  of  the  numerical  boundary  conditions  on  the  edge  of  the  band, 
and  in  this  paper  we  choose  the  boundary  conditions  that  arc  easiest  to  calculate. 

Consider  solving  the  second  order  problem,  equation  (7),  in  Sc.  Since  we  extend  u  off  of  S  in  such  a 
way  that  Vu  •  Vcf)  =  0,  Neumann  conditions  would  be  the  natural  boundary  conditions  to  apply.  However 
for  the  analogous  fourth  order  problem,  Neumann  conditions  are  trickier  to  apply,  especially  on  curved 
domains.  Instead,  for  that  problem,  we  pick  Dirichlet  boundary  conditions,  which  arc  easier  to  implement 
for  the  fourth  order  problem.  Moreover,  Dirichlet  boundary  conditions  assure  us  of  a  symmetric  matrix 
when  we  invert  the  linear  biharmonic  operator,  thus  allowing  us  to  use  the  Jacobian  Conjugate  Gradient 
method.  We  prescribe  u  and  A$u  on  the  boundary,  which  is  more  natural  in  computations  than  prescribing 
u  and  An  on  the  boundary  -  otherwise  we  will  be  required  to  compute  A $u  from  An.  Our  choice  for  A $u  on 
the  boundary  actually  results  from  extending  the  initial  condition  far  enough  past  the  boundary  to  compute 
a  fi  nite  difference  approximation  of  A,n. 

The  numerical  values  of  n  near  the  band  boundary  arc  often  strongly  affected  by  the  boundary  con¬ 
ditions.  We  signifi  cantly  reduce  this  effect  by  regularly  re-extensions  of  n  off  of  S  and  re-initializing  the 
Dirichlet  conditions  based  on  this  re-extension.  By  re-extending  u  up  to  the  boundary,  we  will  have  a  more 
appropriate  boundary  condition  of  u  for  the  next  time-step.  By  re-extending  u  at  least  one  grid  cell  beyond 
the  boundary,  we  can  compute  more  appropriate  values  of  A $u  for  the  next  time  step.  In  the  special  case 
of  the  annular  geometry  of  in  Fig.  1 ,  fi  xing  u  and  A $u  at  the  boundary  is  equivalent  to  fi  xing  u  and  uqq  in 
polar  coordinates. 
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3.7  Visualization 


We  use  MATLAB  to  visualize  the  PDE  dynamics  in  all  of  our  examples.  The  MATLAB  routine  isosurface 
provides  an  easy  means  of  interpolating  values  of  <|)  to  produce  the  surface  S  and  rendering  a  color  map  of 
u  on  S. 


4  Time  Stepping 

A  number  of  choices  arc  available  for  time  integration  of  the  fi  nite  difference  schemes,  each  with  its  own 
advantages  and  disadvantages.  Explicit  schemes  arc  often  used  for  second  order  diffusion  equations  since 
they  arc  easy  to  implement;  to  date  they  arc  the  most  commonly  used  schemes  for  these  surface  problems 
[2,  6].  For  fourth  order  diffusions,  however,  stability  considerations  restrict  explicit  time  steps  to  be  on  the 
order  of  h4.  This  severe  restriction  prohibits  calculating  on  fi  ne  meshes  by  explicit  schemes.  We  thus  turn 
to  implicit  schemes  which  have  no  time  step  restriction  for  linear  fourth  order  diffusion,  though  solving  the 
lineal-  system  may  require  extensive  computation.  Inverting  the  matrices  for  solving  fourth  order  diffusions 
in  three  dimensions  is  computationally  expensive  due  to  the  33  point  stencil  of  the  biharmonic  operator. 
Nonlinear  systems  are  even  more  diffi  cult  to  solve,  as  they  require  repeated  Newton  iterations,  and  time 
steps  may  be  restricted  by  convergence  of  the  Newton  method. 

All  the  above  issues  come  into  play  when  picking  a  method  for  time  stepping.  We  seek  a  method 
that  is  effective  for  strongly  degenerate  PDEs  while  being  easy  to  implement  and  robust  to  changes  in 
both  the  PDE  and  computational  domain.  These  goals  lead  us  to  convexity  splitting  methods,  which  have 
proved  successful  for  the  Cahn-Hilliard  equation  [16],  degenerate  fourth  order  diffusions  including  models 
for  thin  films  [18],  and  for  surface  diffusion  [51],  and  for  second  order  diffusions  on  surfaces  [57],  Our 
convexity  splitting  schemes  require  the  inversion  of  a  linear  system  at  each  time  step.  We  consider  two 
methods  for  inverting  these  systems:  the  conjugate  gradient  method  and  an  Alternating  Direction  Implicit 
(ADI)  method.  Each  method  is  advantageous  for  particular  surfaces,  as  will  be  discussed  below. 

In  this  section,  we  first  describe  convexity  splitting  methods,  then  demonstrate  their  application  to 
lineal'  diffusion  with  Jacobian  Conjugate  Gradient  for  solving  the  linear  systems.  For  now  we  restrict  our 
discussion  to  “simple”  surfaces  with  relatively  small  curvature.  In  section  6  we  turn  to  more  complicated 
surfaces,  and  discuss  how  we  use  ADI  methods  to  invert  the  linear  systems  arising  from  solving  PDEs  on 
these  surfaces. 


4.1  Convexity  Splitting 


Given  a  gradient  fbw  PDE  of  the  form 


ut  —F(u), 


we  split  it  into  two  parts 


ut  =  F\(u)  +F2{u), 


(13) 


where  F\  has  a  strictly  convex  energy  and  /T  has  a  strictly  concave  energy.  Note  this  decomposition  is 
non-unique.  We  update  each  time  step  by  integrating  F\  (u)  implicitly  and  /^(w)  explicitly;  for  example, 
we  might  use  the  scheme 


un+1=un  +  k(Fl(un+l)+F2(un)), 


(14) 


where  w"  denotes  the  value  of  u  after  n  time  steps  of  size  k.  Scheme  (14)  is  unconditionally  stable  in  most 
cases,  and  it  has  an  0{k)  error  [16].  For  gradient  fbws,  this  low  time  stepping  accuracy  is  not  an  issue, 
since  the  objective  is  to  reach  a  local  minimum  of  the  energy.  Although  (14)  may  be  unconditionally 
stable,  larger  time  steps  do  not  necessarily  bring  the  solution  to  this  minimum  any  faster,  rather  there  is 
typically  an  optimal  fi  nite  k  for  obtaining  a  local  minimum  of  the  gradient  fbw’s  energy. 
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This  same  method  may  be  applied  to  problems  that  arc  not  a  gradient  descent;  for  example,  consider 
the  equation 

ii[  =  — V  •  (a  (x)  y  Au). 

Instead  of  computing  this  equation  directly,  we  rewrite  it  as 

u,  —  —CAAu  +  V  •  ((C  —  a  (x))  VA u) .  (15) 

If  C  is  chosen  large  enough,  then  we  can  set 


F\{u)  —  —CAAu 


and 

F2{u)  =  V  •  ((C  —  a  (x))V Au) , 
to  derive  the  unconditionally  stable  scheme 

U  +  ~  U  =  -CAAun+l  +  V  ■  ((C-  a  (x))  VA u”) .  (16) 

K 

The  same  method  has  been  used  for  other  fourth  order  PDEs,  including  surface  diffusion  [51],  the 
Cahn-Hilliard  equation  [16,  54],  and  Hele-Shaw  fbw  [18].  Xu  and  Zhao  used  the  same  idea  to  compute 
second  order  diffusion  on  surfaces  in  [57].  We  describe  exactly  how  C  is  chosen  for  each  of  our  numerical 
examples  in  later  sections.  The  different  equations  we  consider  each  require  different  choices. 

Our  main  challenge  is  to  solve  (16)  quickly  with  a  method  applicable  to  non-rectangular  domains. 
Since  they  computed  in  rectangular  domains,  the  authors  of  [18]  and  [51]  were  able  to  use  Fast  Fourier 
Transform  (FFT)  methods  effectively.  We  experimented  with  similar  methods  for  our  particular  problem, 
and  found  computing  in  a  three  dimensional  box  with  FFT  to  be  slower  than  using  the  methods  we  present 
here  for  computing  in  a  narrow  band  around  the  surface.  Depending  on  the  complexity  of  the  surface 
under  consideration,  we  use  one  of  two  different  methods  for  the  linear  inversion  -  Conjugate  Gradient  or 
Aternate  Direction  Implicit. 

4.2  Example:  Linear  Fourth  Order  Diffusion 

We  fi  rst  demonstrate  convexity  splitting  methods  applied  to  linear  fourth  order  diffusion  on  surfaces, 

ut  =  —A  2su.  (17) 


Assume  |  V(|)|  =  1  and  let  P  denote  the  projection  matrix  of  Section  2,  so  that  (17)  can  be  rewritten  as 

u,  =  -  V-  (PVA su). 

Noting  that  P  =  I  —  N  where  N  is  a  matrix  projecting  onto  the  surface  normal,  V<|),  we  see 

a|m  =  A2u  —  V  ■  (NVA^-m)  .  (18) 

Equation  (18)  also  results  from  setting  C  =  1  in  (15).  At  each  time  step,  we  integrate  the  fi  rst  term  on  the 
right  of  (18)  implicitly  and  the  remaining  term  explicitly: 

un+ 1  _  ,.n 

- - - - =  —A V+1  +  V  •  (NVA sun) .  (19) 

K 
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To  simplify  (19),  we  use  an  approach  suggested  in  [56].  Let  L  =  (/  +  kA2),  subtract  Lun  from  both 
sides  of  (19),  and  group  terms  to  get 

L(un+1  -  un)  =  —kA2u"  +  kV  •  (NVA sun) .  (20) 

We  deb  ne  i-'+l  =  u"+l  —  un  and  notice  the  right  side  of  (20)  is  —kA2u"  to  derive  the  Oik)  scheme, 

Lvn+l  =  -kAjun.  (21) 

When  solving  with  Dirichlet  boundary  conditions,  v  =  0  on  the  boundary,  thus  simplifying  L.  Boundary 
conditions  only  affect  the  explicit  term,  —A 2su'\  which  suggests  fixing  Agu  on  the  boundary,  instead  of 
Am.  Note  that  solving  the  original  form  (19)  requires  using  the  boundary  values  of  both  Am  and  Agu:  the 
boundary  values  of  Am  arc  required  to  defi  ne  L  near  the  boundary  and  Agu  is  needed  to  compute  the  explicit 
term  near  the  boundary.  We  discuss  the  discretizations  of  all  the  spatial  operators  in  Appendix  A.  As  the 
reader  will  notice,  we  use  standard  fi  nite  difference  stencils. 

4.3  Iterative  Solver:  Conjugate  Gradient  Method 

On  domains  with  Dirichlet  boundary  conditions,  the  standard  stencil  for  A2m  is  symmetric,  so  we  can  solve 
(21)  using  the  Conjugate  Gradient  Method.  We  use  the  Conjugate  Gradient  solver  provided  by  ITPACK 
[28].  Unfortunately  a  large  number  of  iterations  may  be  required  at  each  time  step,  especially  for  the 
typically  large  computations  in  three  space  dimensions.  In  cases  where  time  steps  arc  not  limited  by 
stability,  the  large  number  of  iterations  is  not  prohibitive,  and  this  method  is  a  tremendous  improvement 
over  explicit  schemes.  However,  if  the  time  steps  taken  must  be  very  small  (as  we  found  to  be  the  case 
for  surfaces  with  regions  of  high  curvature),  this  method  provides  little,  if  any,  improvement  over  explicit 
methods,  which  require  no  iteration  at  each  time  step. 

5  Numerical  Examples:  Linear  Diffusion  on  Simple  Surfaces 

In  this  section  we  present  the  fi  rst  example,  and  the  simplest  one,  linear  fourth  order  diffusion  on  simple 
surfaces.  Our  basic  computational  framework  follows  the  following  algorithm: 

1.  Store  values  of  the  signed  distance  function  <|)  on  a  Cartesian  grid  in  a  rectangular  domain  containing 
S. 

2.  Mark  grid  points  x/j  (array  entries)  satisfying  |(j(v,/)  |  <  c  for  the  user-defi  ned  bandwidth  c. 

3.  Store  an  initial  value  no  that  satisfi  es  =  /  on  S. 

4.  Evolve  Mo  by  solving  (12)  until  (11)  is  satisfi  ed  in  some  domain  containing  the  band  (far  enough  to 
compute  Aguo  on  the  band’s  boundary). 

5.  Compute  values  of  Aguo  on  the  band’s  boundary. 

6.  Do  while  /"  <  tuax 

(a)  Solve  for  m"+1  in  (21)  using  Jacobi  Conjugate  Gradient  method. 

(b)  Solve  (12)  to  re-extend  values  of  u  off  of  surface  past  band  boundary. 

(c)  Recompute  Agu  on  boundary. 

We  now  present  computational  results  for  (19)  on  both  the  unit  circle  in  M2  and  the  unit  sphere  in  M3.  We 
remind  the  reader  that  discretizations  of  all  spatial  derivatives  arc  described  in  Appendix  A. 
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5.1  Diffusion  on  the  Unit  Circle 

We  compare  our  computations  with  an  exact  solution  of  (17)  on  the  unit  circle.  In  this  case,  solving  (17) 
is  exactly  the  same  as  solving 


Ut  —  U.SSSS  (22) 

«M)  =  f(s)  (23) 

on  [0. 2n]  with  periodic  boundary  conditions.  Equation  (22)  can  be  solved  exactly  using  separation  of 
variables.  We  choose  f(s)  =  cos(4.v).  giving  the  exact  solution 

u(x,t)  =  cos(45j<?-64'. 

We  use  conjugate  gradient  method  to  invert  the  implicit  term.  We  compare  our  results  with  the  exact  value 
of  u  at  times  t  =  0.001  and  t  =  0.01.  Although  these  times  arc  small,  the  half  life  of  (17)  for  this  initial  data 
is  approximately  t  =  0.01.  We  defi  ne  the  level  set  function  on  the  box  [—2,2]  x  [—2,2]  using  grid  sizes  of 
h  =  where  M  =  40, 80, 160,  and  320.  We  compute  in  a  band  of  width  3/8  in  these  calculations.  Since 
our  scheme  is  0(k),  we  use  a  time-step  of  k  =  -y^/r.  We  must  choose  a  small  constant  of  proportionality  to 
get  adequate  results.  We  expect  this  is  due  to  the  high  degeneracy  of  the  equation.  We  must  choose  a  small 
constant  of  proportionality  to  get  adequate  results.  We  expect  this  is  due  to  the  high  degeneracy  of  the 
equation.  The  error  in  the  long-time  calculation  is  of  course  improved  by  careful  re-extension,  although 
we  report  here  the  results  without  it  in  order  to  better  identify  the  contributions  of  the  different  components 
of  our  scheme.  Note  that  for  M  =  320,  over  300  time  steps  arc  necessary  to  reach  t  =  0.01.  Without  re¬ 
extension  of  data  off  of  the  surface,  the  Dirichlet  boundary  conditions  increase  the  solution’s  error.  Figure 
2  demonstrates  an  error  of  ()(h2)  in  the  short-time  calculation  and  O(h)  in  the  long-time  calculation.  This 
decrease  in  accuracy  is  caused  by  the  choice  of  Dirichlet  Boundary  conditions.  We  found  that  we  maintain 
0(h2)  accuracy  by  re-extending  the  initial  data  after  every  five  time  steps.  To  calculate  the  error,  we 
interpolated  to  fi  nd  the  values  of  u  on  the  circle  and  then  used  the  t°  norm  of  the  difference  of  these  values 
from  the  exact  solution. 

5.2  Diffusion  on  the  Unit  Sphere 

Figure  3  shows  linear  fourth  order  diffusion  on  the  unit  sphere.  Implementation  requires  only  slight 
changes  from  the  example  of  a  circle.  We  only  need  to  change  the  level  set  function  and  use  fi  nite  dif¬ 
ferences  in  three  dimensions.  This  easy  adaptability  to  higher  dimensions  is  typical  of  level  set  methods. 
In  addition,  we  may  easily  compute  on  surfaces  such  as  ellipsoids  or  tori  by  changing  only  the  level  set 
function. 

In  our  example,  the  level  set  function  (|)  isdefi  ned  on  [—2,2]  x  [—2,2]  x  [—2,2]  with  100  x  100  x  100 
grid  points.  The  initial  condition  is  defi  ned  using  polar  coordinates  by  rq(p. 0,  (3)  =  sin (30)  sin (7(3).  We 
take  time  steps  k  =  5 h2  where  h  is  the  grid  cell  width.  Figure  3  displays  the  solution  at  t  =  0.  t  =  0.2, 
and  t  —  0.9.  We  use  a  bandwidth  of  ten  grid  cells  (fr  ve  cells  off  the  surface  in  each  direction),  and  re¬ 
extend  the  initial  data  every  four  time  steps  by  solving  (12)  to  steady  state.  See  Appendix  A  for  the  spatial 
discretizations  of  the  operators  in  (21). 

6  ADI  and  Complex  Geometry  Surfaces 

In  Section  5,  we  solved  linear  diffusion  on  two  simple  geometries:  the  unit  circle,  and  the  unit  sphere.  For 
these  examples.  Conjugate  Gradient  Method  provided  an  effi  cient  means  of  inverting  the  linear-  systems 
in  our  implicit  schemes.  Unfortunately,  as  we  discovered  experimentally,  the  Conjugate  Gradient  Method 
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Error  at  t  =  0.001 


Error  at  t  =  0.01 


Figure  2:  Error  plot  for  linear  fourth  order  diffusion  on  the  unit  circle,  without  re-extension. 


often  requires  a  large  number  of  iterations.  This  becomes  prohibitive  when  forced  to  take  small  time 
steps,  as  we  found  to  be  the  case  for  solving  PDEs  on  surfaces  with  regions  of  high  curvature  (see  Section 
3.1).  For  these  surfaces,  we  turn  to  Alternate  Direction  Implicit  (ADI)  schemes,  which  do  not  require 
repeated  iterations.  ADI  schemes  involve  the  inversion  of  only  banded  matrices,  which  require  only  0(M) 
operations  for  an  M-by-M  matrix.  ADI  schemes  have  been  used  for  both  second  and  fourth  order  PDEs 
[12,  13,  14,  56,  58].  For  example,  assume  we  arc  solving  the  heat  equation  on  a  rectangle  with  an  implicit 
scheme: 


,n+ 1 


—  •/'+' 

*xx 


,.n+ 1 

*yy  ■ 


(24) 


The  following  is  an  Oik)  approximation  of  (24): 


(I  -kd2x)(I-kd*)un+1  =un. 


(25) 


Unlike  (24),  equation  (25)  only  requires  the  inversion  of  tridiagonal  matrices. 

Lineal-  fourth  order  diffusion, 

iii  =  —A  2u,  (26) 

is  more  diffi  cult  to  implement  with  ADI,  as  A2  includes  a  cross-term,  23232.  In  [56],  Witelski  and  Bowen 
suggest  an  ADI  scheme  in  which  the  mixed  derivative  term  is  computed  explicitly.  They  showed  that 

{I  +  kd*){I  +  kd*)un+1  =  (I  — 2k3232)un  (27) 

is  an  unconditionally  stable  0(k)  scheme  for  (26).  Higher  order  accurate  ADI  schemes  following  from  the 
same  idea  are  discussed  in  [56]. 

Unfortunately,  applying  an  ADI  scheme  directly  to  the  degenerate  diffusions  we  consider  here  yields  a 
scheme  with  a  stability  restriction  that  is  no  better  than  for  explicit  schemes.  The  interested  reader  will  fi  nd 
a  discussion  of  this  in  Appendix  B.  In  order  to  improve  upon  the  stability  restrictions  of  explicit  schemes, 
we  combine  ADI  with  convexity  splitting. 
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Figure  3:  Linear  fourth  order  diffusion  on  the  sphere  after  0,  4,  and  18  time  steps.  Red  denotes  u  =  1  and 
blue  denotes  u  —  —  1 . 


16 


6.1  Example:  Linear  Fourth  Order  Diffusion 

As  a  fi  rst  example,  consider  solving  (19)  with  an  ADI  method.  Rewrite  (19)  as 

(/  +  kA2)un+1  =  un  +  JfcV  •  (NVA5n") .  (28) 

Letting  Lx  —  (7  +  kdx)  and  Ly  =  (/  +  kc)4).  we  introduce  O(k)  errors  by  modifying  (28)  to 

LxLyun+ 1  =  (/  -  2 kd2d2)un  +  kV  •  (NVA sun ) .  (29) 

Now  subtract  LxLyun  from  both  sides  and  deb  ne  ri+1  —  v"  to  derive  the  compact  scheme 

LxLyvn+ '  =  -A  2sun.  (30) 

Unlike  Locally  One-Dimensional  (LOD)  schemes,  this  method  easily  extends  to  three  dimensions, 
producing  an  unconditionally  stable  scheme, 

LxLyLzv"+l  =  -A  2su\  (31) 

with  Lz  —  (/  +  kd4 ) .  This  easy  application  to  three  dimensions  is  obviously  important  for  the  problems  we 
consider.  However  ADI  schemes  have  a  disadvantage  in  non-rcctangular  domains  such  as  the  narrow  band 
(|<K*)I  —  c)  we  consider.  This  is  best  seen  by  rewriting  (27)  as 

W  =  LyUn+1 

Lxw  =  {I-2kd2xd2)un.  (32) 

At  the  half  time  step  we  solve  for  w.  which  does  not  have  the  same  boundary  conditions  as  u.  Solving 
for  these  boundary  conditions  on  non-rcctangular  domains  is  a  tricky  problem  that  becomes  more  diffi  cult 
for  complicated  domains  -  domains  like  the  bands  we  wish  to  compute  in.  Instead  we  approximate  the 
boundary  condition  by  fixing  w  —  u  on  the  boundary,  as  suggested  in  [58].  Unfortunately  this  method 
is  not  unconditionally  stable  as  we  discovered  in  our  numerical  experiments.  Even  with  this  restriction, 
however,  we  found  this  ADI  method  to  be  more  effi  cient  than  Conjugate  Gradient  Method  for  solving 
PDEs  on  surfaces  with  complicated  geometries. 

6.2  Numerical  Example:  Linear  Fourth  Order  Diffusion  on  the  Stanford  Bunny 

The  bunny  data  for  this  example  is  taken  from  [30].  Due  to  the  existence  of  regions  of  high  curvature  on 
the  bunny,  we  fi  rst  run  the  standard  heat  equation  in  if  on  the  bunny’s  level  set  function  for  a  very  short 
time  period.  This  smoothes  any  singularities  in  the  level  set  function  for  the  bunny  while  maintaining  its 
essential  structure,  as  Figure  4  shows.  We  found  this  surface  smoothing  to  be  unnecessary  for  solving 
second  order  equations  on  the  bunny  -  it  is  required  for  higher  order  equations  at  the  grid  resolution  with 
which  we  work.  After  smoothing  <|),  we  can  either  use  this  new  (|)  as  is,  or  we  may  keep  only  its  zero  level 
set  then  reinitialize  (|>  so  that  it  is  again  a  distance  function. 

We  follow  the  algorithm  in  Section  5  with  only  two  modifi  cations.  As  just  discussed,  there  is  an 
added  preprocessing  step  where  the  level  set  function  for  the  bunny  is  smoothed,  and  instead  of  solving 
our  numerical  scheme  with  the  Conjugate  Gradient  method  at  each  time  step,  we  use  the  ADI  scheme 
(31).  Figure  4  shows  fourth  order  linear  diffusion  on  the  surface  of  the  bunny.  We  use  a  cubic  grid 
whh  cell  width  li  =  |.  There  arc  159  x  161  x  129  grid  cells.  We  use  a  bandwidth  of  about  10  grid  cells. 
The  time  step  is  k  =  ^L_/z3.  When  computing  with  an  explicit  scheme,  we  needed  k  =  ^y/;4,  so  there 
is  still  signifi  cant  improvement  over  the  explicit  method.  We  defi  ne  the  initial  condition  in  all  of  R  by 
uo{x,y,z)  =  ^{cos{ |x)  +  .sin(ly)).  We  then  fi  x  u  on  the  bunny  and  use  the  extension  procedure  to  change 
the  values  of  u  off  of  the  surface  so  that  the  initial  condition  satisfies  V/q  •  V()  =  0.  We  re-extend  every 
four  time  steps. 
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t  —  0.2 

Figure  4:  Linear  fourth  order  diffusion.  Each  row  shows  two  views  of  the  surface  at  the  same  point  in 
time. 
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7  Cahn-Hilliard  Equation 


We  now  consider  a  classical  phase-fi  eld  model  for  spinodal  decomposition  of  alloys,  the  Cahn-Hilliard 
equation  [41].  Such  models  describe  coarsening  dynamics  such  as  the  phase  separation  following  a  quench 
from  a  disordered  to  an  ordered  phase.  Computer  simulations  play  an  important  role  in  the  characterization 
of  late-stage  coarsening  processes.  Recent  efforts  have  focused  on  developing  numerical  methods  for 
the  Cahn-Hilliard  equation  in  Euclidean  geometries  using  fi  nite  element  methods  [4]  and  psuedospectral 
methods  [54].  Here  we  develop  a  numerical  method  for  solving  the  Cahn-Hilliard  equation  on  a  general 
surface,  using  the  implicit  representation  approach. 

We  solve 

ut  —  A s  (— £2As  u  +  u3  —  u) ,  (33) 

for  which  u  —  1  and  u  —  —  1  are  both  stable  steady  states.  In  our  examples  we  use  an  initial  condition  of 
u  =  0  plus  a  very  small  zero  mean  perturbation.  The  solution  u  quickly  separates  the  surface  S  into  two 
regions  S+  and  S-  where  u  takes  on  values  of  1  and  —  1  respectively.  The  remaining  points  of  S  lie  on  the 
interface  of  width  0(e)  between  these  two  regions.  In  later  stages,  u  undergoes  spinodal  decomposition; 
S. |_  and  S-  change  shape  so  that  the  length  of  the  interface  between  the  two  regions  decreases  while 
maintaining  the  area  of  each  region  (and  the  mean  value  of  u).  This  coarsening  of  S'-),  and  S_  slows  with 
time. 

In  free  space,  the  Cahn-Hilliard  equation  is  a  diffuse  interface  model  for  Mullins-Sekerka  dynamics 
[40,  44],  and  a  version  of  the  Cahn-Hilliard  equation  with  degenerate  mobility  is  a  diffuse  interface  model 
for  Hele-Shaw  [18],  For  fbws  on  surfaces,  we  might  expect  to  be  able  to  use  these  models  as  diffuse 
interface  models  of  curve  evolution  on  complicated  geometries.  This  idea  is  further  suggested  by  the 
simulations  shown  in  fi  gures  5  and  6. 

Equation  (33)  is  slightly  more  complicated  than  (17)  due  to  the  nonlinear  second  order  term.  We  pick 
Fx(u)  =  — e2  A2  u  in  (13)  to  derive  the  scheme 

- — —  —  =  — e2A 2un+l  +e2V  •  (NVAsun)  +  AS  ((rr")3  -  n"j  ,  (34) 

which  has  a  stability  requirement  of  k  —  0(h2).  We  defi  ne  1  =  un+l  —  un  and  L  —  I  +  ke2A2  to  simplify 
(34)  as  was  done  in  Section  4.2: 


Lvn+1  =  As(-e2A sun  +  (un)3  -  un).  (35) 

Note  that  our  method  does  not  split  the  energy  for  (33)  into  convex  and  concave  parts.  Doing  so  requires 
inverting  an  operator  involving  A<)m"+i  ,  and  we  found  this  to  be  prohibitively  slow  in  practice. 

We  fi  rst  solve  (33)  on  a  torus.  The  surface  does  not  seem  to  affect  the  stability  of  (35),  so  we  take  large 
time  steps  (relative  to  explicit  schemes)  and  solve  the  linear  systems  with  Jacobi  Conjugate  Gradient,  as 
discussed  in  Section  4.3.  We  choose  £  =  2 h  and  fi  x  the  initial  condition  to  he  /q  =  0 + q ,  where  r|  is  a  small 
random  perturbation  across  the  surface.  Our  choice  of  £  is  similar  to  that  used  in  [54]  for  Cahn-Hilliard 
with  constant  mobility  and  in  [18]  for  Cahn-Hilliard  with  degenerate  mobility.  Our  results  have  transition 
layers  approximately  four  to  six  grid  cells  thick.  We  compute  on  a  120  x  120  x  120  grid  with  cell  width 
h  —  y  The  time  step  taken  is  k  —  j^h2.  We  observe  the  spinodal  decomposition  expected;  see  Figure  5. 

Figure  6  shows  spinodal  decomposition  in  the  Cahn-Hilliard  equation  on  the  smoothed  bunny.  We 
use  ADI  to  compute  (35)  as  in  Section  6.  We  arc  constrained  to  take  the  same  time  step  as  in  the  linear 
problem,  k  =  0(h3).  Again  £  =  2 h  and  u$  =  0  +r|,  where  r|  is  a  small  zero  mean  perturbation. 
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t  =  0.45 


t  =  3.5 


Figure  5:  Spinodal  decomposition  in  the  Cahn-Hilliard  equation.  Initial  condition  is  u  —  0  (slightly  per¬ 
turbed),  red  denotes  u  =  1  and  blue  denotes  u  =  —  1.  Each  row  shows  two  views  of  the  surface  at  the  same 
point  in  time. 
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t  =  950 


Figure  6:  Spinodal  decomposition  in  the  Cahn-Hilliard  equation.  Initial  condition  is  u  —  0  (slightly  per¬ 
turbed),  red  denotes  u  =  1  and  blue  denotes  u  =  —  1.  Each  row  shows  two  views  of  the  surface  at  the  same 
point  in  time. 
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8  Thin  Film  Fluid  Flow 


An  interesting  class  of  problems  described  by  nonlinear  fourth  order  diffusion  equations  arc  thin  fi  lm 
fbws  involving  a  layer  of  viscous  liquid  on  a  solid  surface.  Surface  tension  forces  lead  to  the  fourth  order 
motion,  due  to  curvature  affects  on  the  air-liquid  interface  of  the  fi  lm  [7,  34,  42].  In  the  case  of  such  fbws 
on  curved  surfaces,  the  underlying  curvature  of  the  substrate  plays  a  role  in  the  motion  of  the  fi  lm,  with 
fbw  out  of  regions  of  high  curvature.  Recently  there  has  been  some  interest  in  the  fLiids  community  in 
deriving  equations  of  motion  for  such  fbws  on  general  surfaces.  The  most  comprehensive  paper  to  date  on 
this  problem  is  that  of  Roy,  Roberts  and  Simpson  [47],  building  on  previous  work  of  Schwartz  and  Weidner 
[48].  Several  authors  have  developed  numerical  methods  for  thin  films  on  surfaces  using  coordinates  on 
the  surface.  This  work  includes  fbws  on  cylinders  [24,  48,  55]  and  more  general  surfaces  [37],  as  well  as 
icing  of  airplane  wings  [39,  35,  36,  38].  We  develop  a  numerical  method  for  thin  fi  1ms  on  general  surfaces 
using  the  implicit  representation  methodology. 

Letting  u  denote  the  fi  lm  thickness,  Roy,  Roberts,  and  Simpson  derived  the  following  for  the  dimen¬ 
sionless  fbw  of  a  thin  fi  lm  in  the  absence  of  gravity  [47] : 


dA 

dt 


u2tysk  -  - u 4  (kI  -  K)  •  V5k 


(36) 


K  is  the  curvature  tensor  for  the  curved  substrate.  The  quantities  k  \  and  k2  arc  the  principle  curvatures  of 
the  substrate,  and  k  =  k\  +k2  is  twice  the  mean  curvature,  k  is  curvature  of  the  free  surface,  which  we 
approximate  as 

k  =  K  +  (kj  +  +  A  $u. 

The  variable  £  is  the  amount  of  find  above  a  surface  patch,  which  we  approximate  by  C,  =  u  —  \ku2  + 

Using  the  framework  of  our  level  set  function  (|),  we  can  easily  compute  K  and  K,  as  well  as  the 
quantities  k\.k2.  and  k\  +  k\.  We  note  that  K  is  the  Jacobian  of  the  Gauss  Map  of  the  substrate.  In  terms 
of  the  implicit  representation,  the  Gauss  Map  is  The  trace  of  K  gives  — k,  and  the  determinant  of 
A  =  K  +  gives  the  Gauss  curvature,  k{k2 .  We  fi  nally  note  that 

k\  +  kl  =  K2  -2k\k2. 


This  model  differs  from  both  linear  diffusion  and  the  Cahn-Hilliard  equation  in  that  the  highest  order 
term  of  (36)  is  nonlinear.  Our  convexity  splitting  must  take  this  nonlinearity  into  account.  To  derive  the 
evolution  scheme,  we  carry  out  the  differentiation  on  the  left  of  (36)  to  get  an  evolution  equation  for  u  : 


dt  3£'  5 


«%V5k--M4(Kl  — K)-Vsk 


where  C,1  —  1  —  Kit  +  K|  K?/r.  Following  our  earlier  examples,  we  pick  an  appropriate  C  and  evolve 


un+l  -  un 


=  —C"AAu"+ 1  +CnAAun  -F{un), 


(37) 


(38) 


where  F(un)  denotes  the  right  hand  side  of  (37).  The  dependence  of  C  on  the  time  step  relicts  the  depen¬ 
dence  of  the  highest  order  term  of  (37)  on  un.  Noting  that  the  highest  order  term  of  F(u")  is  j^ir^A^u11 , 

we  see  that  choosing  Cn  <  ^ u2^  might  cause  the  method  to  be  unstable.  On  the  other  hand,  picking  C  too 
large  signifi  candy  slows  the  dynamics.  So  again  using  A  to  denote  the  computational  band  <  c,  we 

pick 


Cn  =  max 


*,j,kesc  3(1  -  K u'\j  k  +  k1k2(«'!j/) 


«m)2Cuu 


(39) 
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8.1  Example:  Thin  Film  on  an  Ellipse 


The  curvature  dependent  terms  in  (36)  refbct  fluid  motion  driven  by  curvature  of  the  surface.  Fluid  builds 
up  in  regions  of  high  negative  curvature  while  leaving  areas  of  high  positive  curvature.  We  consider  an 
example  discussed  in  [47]  and  study  the  fbw  of  an  initially  constant  layer  of  fi  lm  on  an  ellipse.  The  authors 
of  [47]  discussed  the  plausible  motion  of  the  Ibid,  however  they  did  not  compute  (36)  for  this  example. 
We  provide  numerical  evidence  supporting  their  claims  for  the  fluid  dynamics. 

We  place  a  layer  of  fbid  with  constant  thickness  u  =  0.02  on  the  outside  of  the  ellipse 

x2  +  3y2  =  2.  (40) 

We  computed  on  [—2,2]  x  [—2,2]  with  200  x  200  grid  points.  We  used  a  constant  time  step  k  —  10/? 2 
where  h  denotes  the  grid  cell  width.  We  used  the  algorithm  in  Section  5  with  only  one  change:  C  is 
updated  dynamically  as  described  above.  As  seen  in  Figure  7,  the  Ibid  moves  away  from  the  ellipse’s 
regions  of  high  curvature  near  y  =  0. 

Now  consider  placing  the  same  constant  layer  of  Arid  on  the  inside  of  the  ellipse  (40).  Once  the  im¬ 
plicit  function  (|>  for  the  above  problem  is  defi  ned,  we  easily  adjust  our  simulation  to  this  case  by  mapping 
(|)  — >  — 4>.  On  the  inside  of  the  ellipse,  the  fi  lm  builds  up  in  the  areas  of  high  curvature.  See  Figure  8. 


8.2  Sphere  with  Gravity:  Fingering  Instability 

The  following  model,  derived  in  [47],  includes  the  effects  of  gravity 

1 


*  =  -V 

dt  3  S 

1 


?r i^V^K  —  Tm4  (kI  —  K)  •  Vk 


--BoVs- 


iSgs  -  it  (  Kl  +  )  •  gs  +gnu3Vsu 


(41) 


The  Bond  number  Bo  quantifi  es  the  relative  strength  of  gravity  to  surface  tension,  and  it  is  given  by 


Bo  =  p  gH-a, 


where  p  is  the  fbid  density,  g  is  acceleration  due  to  gravity,  a  is  the  surface  tension,  and  H  is  the  charac¬ 
teristic  thickness  of  the  fi  lm.  The  vector  gs  is  the  component  of  gravity  tangent  to  the  surface  S',  and  gn 
is  the  magnitude  of  the  component  of  gravity  normal  to  S.  In  our  example,  g  =  —  (0,0, 1),  so  gs  =  Py  and 

gn  =  II -gs\ ■ 

We  place  a  ring  of  fbid  with  sinusoidally  varying  height  near  the  top  of  the  sphere.  The  remainder 
of  the  sphere  has  a  precursor  layer  of  thickness  u  —  5  x  10-3.  Gravity  drives  the  fbid  to  the  bottom  of 
the  sphere  while  the  non-constant  thickness  causes  a  fi  ngering  effect.  We  used  up-winding  to  compute  the 
gravity  terms  in  (41).  We  computed  in  [—2,2]  x  [—2,2]  x  [—2,2]  on  a  128  x  128  x  128  grid.  We  used  a 
constant  time  step  k  =  \lr  where  h  denotes  the  grid  spacing.  We  compute  in  a  band  that  is  about  ten  grid 
cells  wide  and  re-extend  data  off  of  the  surface  after  every  fi  ve  time  steps.  The  Bond  number  is  Bo  =  100  in 
this  simulation.  The  results  are  given  in  fi  gures  9,  10,  and  11.  Note  that  the  fbid  drips  down  the  sphere  and 
eventually  begins  to  collect  at  the  bottom  of  the  sphere.  Equation  (41)  is  derived  from  lubrication  theory 
and  thus  does  not  capture  later  time  dynamics  in  which  fbid  might  drop  off  the  bottom  of  the  sphere. 
Although  the  computation  is  performed  on  a  spherical  surface,  the  method  is  very  general  and  can  easily 
be  applied  to  surfaces  that  lack  the  symmetry  of  the  sphere. 
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Film  Outside  of  Ellipse 


angle  of  position  with  x-axis 


Figure  7:  Thin  fi  lm  driven  by  curvature  on  the  outside  of  an  ellipse. 
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Film  Inside  Ellipse 


angle  of  position  with  x-axis 


Figure  8:  Thin  fi  lm  driven  by  curvature  on  the  inside  of  an  ellipse. 
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t  =  0.0 


Figure  9:  Equation  (41)  solved  on  the  sphere.  Each  row  gives  a  side  and  bottom  view  of  the  sphere  at  the 
same  point  in  time.  The  sphere  is  covered  with  a  precursor  layer  of  thickness  of  u  —  5  x  10-3.  The  color 
scheme  is  used  to  display  fi  lm  thickness. 
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Figure  10:  Equation  (41)  solved  on  the  sphere.  Each  row  gives  a  side  and  bottom  view  of  the  sphere  at  the 
same  point  in  time.  The  sphere  is  covered  with  a  precursor  layer  of  thickness  of  u  —  5  x  10-3.  The  color 
scheme  is  used  to  display  fi  lm  thickness. 
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t  =  16.0 


0.01  0.02  0.03  0.04  0,05  0.05  0.07 


Figure  11:  Later  dynamics  of  equation  (41)  solved  on  the  sphere.  Each  row  gives  a  side  and  bottom  view 
of  the  sphere  at  the  same  point  in  time.  The  sphere  is  covered  with  a  precursor  layer  of  thickness  of 
u  =  5  x  10-3.  The  color  scheme  is  used  to  display  fi  lm  thickness. 
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9  Conclusions 


We  have  presented  a  method  for  solving  fourth  order  PDEs  on  surfaces  of  arbitrary  geometry  with  fi  nite 
difference  schemes  on  a  Cartesian  mesh.  With  our  methods,  the  same  code  can  be  easily  applied  to  many 
different  surfaces  by  changing  only  the  function  implicitly  defining  the  surface.  In  computing  higher 
order  PDEs  on  surfaces,  we  derive  equations  that  arc  solved  in  an  augmented  band  around  the  surface  in 
Euclidean  space.  As  PDEs  on  WN .  these  equations  arc  unique  in  the  severity  of  their  degeneracy  -  we 
know  of  no  other  fourth  order  diffusions  being  studied  where  all  diffusion  locally  is  in  one  coordinate 
direction  only. 

The  higher  order  equations  introduce  a  number  of  challenges  that  arc  of  little  consequence  for  fi  rst  or 
second  order  PDEs.  Stability  restrictions  for  time  stepping  fourth  order  equations  require  using  implicit 
schemes,  unlike  most  previous  work  on  solving  PDEs  on  implicit  surfaces  [2,  6,  32],  We  derived  semi- 
implicit  schemes  using  convexity  splitting  ideas  explored  in  [16,  18,  51]  and  presented  a  new  means  of 
combining  convexity  splitting  schemes  with  ADI  methods.  Compared  to  lower  order  equations,  fourth 
order  PDEs  require  more  careful  extension  and  re-extension  of  data  off  of  the  surface,  and  they  have  more 
complicated  boundary  conditions.  We  discussed  each  of  these  issues  in  detail  and  applied  our  methods  to 
lineal-  fourth  order  diffusion,  the  Cahn-Hilliard  equation,  and  a  recently  derived  model  for  surface  tension 
driven  fbws  on  curved  substrates. 

Our  work  is  only  a  fi  rst  step.  Although  our  schemes  are  faster  than  explicit  schemes,  there  is  room 
for  improvement.  It  remains  a  diffi  cult  problem  to  improve  the  time  step  restriction  for  fbws  on  complex 
surfaces  with  high  curvature,  because  the  inherent  geometry  is  embedded  in  the  projection  operator. 

It  is  an  interesting  problem  to  develop  schemes  with  higher  numerical  accuracy,  in  particular  for  the 
time  step,  as  this  might  allow  for  computations  with  a  larger  time  step.  If  the  geometry  introduces  stiffness 
into  the  dynamics,  these  terms  could  also  be  dealt  with  using  splitting  methods.  Local  mesh  refi  nement, 
while  still  taking  advantage  of  the  Cartesian  geometry,  might  also  help  address  this  issue. 

It  is  our  hope  that  the  techniques  developed  in  this  paper  will  be  useful  for  scientists  interested  in 
computing  higher  order  PDEs  on  complicated  geometries.  Thin  fi  lm  fbw  is  a  particularly  relevant  problem 
where  there  is  interest  in  modeling  the  liquid  lining  of  the  lungs,  icing  of  airplane  winds,  and  numerous 
other  problems  in  which  the  underlying  surface  geometry  is  complex. 

A  Spatial  Discretization 

We  now  present  the  specifi  c  discretizations  for  the  gradient,  Laplacian,  and  Laplace-Beltrami  operators. 
We  compute  everything  on  a  Cartesian  grid  in  M;V .  Assume  the  domain  is  given  by  0  <  x/  <  L  for  1  <  i  <  N. 
We  use  a  cubic  grid  with  Ax,-  =  h  for  each  /.  Let  (0  =  (or,)  be  a  vector  of  integers  0  <  or,-  <  j  and  denote 
the  coordinate  directions  by  e,,  so  that  each  grid  point  is  given  by 

*oo  =  Yj  /iC0'e<- 

0<i<N 

Our  discretizations  follow  the  example  of  [6],  We  calculate  V([)  with  centered  differences, 

V<K*co)  =  +  hei)  ~  <K*g>  -  (42) 

In 

but  use  forward  differences  to  compute  gradients  for  everything  else: 

Vw(x oo )  =  j{u(xa  +  hei)  -u(xm))j.  (43) 

n 
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Let  v(xM)  be  a  vector  whose  /th  component  is  v'(xM).  We  compute  the  divergence  of  v  using  backward 
differences: 

V-v(xo,)  =  y-  £  v'(*to)  -v'(xo)-/;e,).  (44) 

,l  0 <i<N 

The  Laplacian  of  a  scalar  quantity  u  is  computed  by  applying  (44)  to  (43),  and  the  biharmonic  is  computed 
by  repeated  application  of  the  Laplacian. 

Projected  derivatives  arc  computed  using  appropriate  averaging  of  the  projection  matrix.  Given  a 
projection  matrix  P  =  (piJ).  we  compute 


PV u(x to)  =  (  Y,  P u*j  (x“) ) 
i<;<w 


where 

ftuXj{xw)  =  —  {pli{xm  +  hej)  +plj(xa))(u(xta+hej)  -  u(x, «,))- 
Assuming  that  cf)  is  a  signed  distance  function,  by  applying  (44)  to  (45)  we  compute 


Asu(x. a,)  =  V-  (PVw(xa))) 


I 

1  <i,j<N 


jjuXj  (xro)  -  PuXj  (xco  -  h&i) 
h 


(45) 


(46) 


We  repeat  (46)  to  calculate  A$A $u. 

B  Degeneracy  and  ADI 

To  demonstrate  the  strong  degeneracy  of  the  fourth  order  equations  considered  in  this  paper,  we  apply 
standard  ADI  schemes  to  both  (7)  and  (9)  without  the  assistance  of  the  operator  splitting  discussed  in 
Section  4.1.  Our  calculations  suggest  the  need  for  an  indirect  method  (such  as  convexity  splitting),  while 
also  showing  that  the  second  order  problem  is  far  less  delicate  than  its  higher  order  counterpart. 

We  examine  the  special  case  where  S  is  the  line  in  M2  making  an  angle  of  0  with  the  x-axis.  We  fi  rst 
discuss  the  ADI  method  applied  to  second  order  diffusion  on  the  line,  then  continue  with  a  discussion  of 
ADI  for  intrinsic  fourth  order  diffusion  on  the  same  line.  We  show  that  the  stability  of  the  fourth  order 
diffusion  scheme  depends  on  0,  which  hardly  affects  the  second  order  diffusion  scheme.  For  simplicity  we 
discuss  only  the  two  dimensional  problem,  but  note  that  similar  results  hold  for  three  dimensions. 

B.l  ADI  without  Convexity  Splitting 

To  apply  ADI  to  equations  (7)  and  (9)  without  convexity  splitting,  we  follow  the  example  of  [56]  and 
compute  terms  with  mixed  derivatives  explicitly  and  all  remaining  terms  implicitly.  We  fi  rst  examine 
equation  (7),  leaving  equation  (9)  for  Section  B.2.  Letting  S  be  the  line  making  angle  0  with  the  x-axis, 
our  distance  function  is 

(|)  =  —  sin  0  x  +  cos  0  y. 

and  the  projection  matrix  is 

_  cos20  cos  0  sin  0 

cos  0  sin  0  sin20 

For  this  example,  (5)  simplifi  es  to 

ut  =  cos20  nxr  +  2cos0sin0  +  sin20  uyy.  (48) 
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Defi  ne 


and 


Dx  -  cos2  0  c)xx. 
Dy  —  sin2  0  dyy, 
Lx  —  I  kDx, 


Ly  =  I  —  kDy. 


Note  that  LxLyiin  =  (/  —  kDx  —  kDy  +  k2DxDy)  so  ignoring  the  term  of  0(k2) ,  we  see  that  LxLy  approx¬ 
imates  (48)  without  its  mixed  partials.  We  now  introduce  O(k)  errors  by  computing  the  mixed  partials 
explicitly.  This  can  be  done  more  accurately  so  that  these  errors  arc  reduced  -  see  [56].  So  we  have 

LxLy  u“-'  =  — 2k  cos  Osin  0  unxy. 

We  know  subtract  LxLyu"  from  both  sides  of  the  above  equation  and  let  v',+1  =  u”+l  —  u"  to  get 

LxLyvn+ 1  =  kV  •  (P Vw" )  =  k  (cos2  0  unxx  +  2  cos  0  sin  0  +  sin2  0  unyy ) .  (49) 

We  use  Von  Neumann  stability  analysis  to  study  scheme  (49)  for  this  example.  Solutions  of  the  con¬ 
tinuous  in  space  and  discrete  in  time  solutions  of  (49)  can  be  written  as 


un(x,y)  =  a(a)"exp/(aix  +  a2y), 


(50) 


where  the  superscript  n  on  the  left  corresponds  to  the  n-th  time  step,  while  the  superscript  on  the  right 
corresponds  to  a  power  of  o.  Following  the  example  of  [56],  we  rescale  the  wave  number  as  a  =  ia  and 
substitute  (50)  into  (49)  to  get 

,  (cosOoci  +  sin0ccT)2 

c(a)  -  1  = - - — 2 - ^ ^ - V - - — z — r. 

1  +  cos-  0oc|  +  sin"  Oa^  +  sin-  0  cos-  OocjoCt 

We  see  (49)  is  unconditional  stable,  since 

(cosOoci  +sin0(X2)2 


1  +  cos2  0ay  +  sin2  +  sin2  0cos2  0a2 


<  2 


(51) 


is  satisfi  ed  for  all  a  and  0.  Though  simpler  than  (7)  on  general  surfaces,  this  example  shows  that  (49)  is 
stable  even  when  the  diffusion  is  not  in  the  direction  of  the  x  or  y  axes.  We  next  show  that  the  fourth  order 
problem  is  very  different,  as  it  is  too  degenerate  for  this  type  of  ADI  scheme. 


B.2  Fourth  Order  Problem 

We  now  discretize  (9)  as  in  [56].  Defi  ne 


and 


Dx  —  cos  0  ()  vvvv . 

Dy  -  sin  0  rlyyyy , 

Lx  =  I  +  kDx, 

Ly  =  I  +  kDy. 
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Once  again  letting  v‘ 


—  un.  the  scheme 


"+t  =  un+l 


LxLyvn+]  =  -kV  •  (PV  •  (PVn”)) 


(52) 


is  an  0(k)  approximation  of  (9). 

We  repeat  the  Von  Neumann  analysis  used  for  (49),  this  time  rescaling  a  =  p-a  to  get 


a(  a)  - 1  =  - 


(cos0ai  +  sin0ct2)4 

1  +  cos4  0a4  +  sin4  0a4  +  sin4  0  cos4  0a{  a4 


Stability  of  (52)  thus  requires 


_ (cos0oti  +  sin0q2)4 _  <  ^ 

1  +  cos4  0a4  +  sin4  0a4  +  sin4  0  cos4  0a4  a4  — 

Note  that  (53)  is  clearly  satisfied  when  sin0  =  0  or  cos0  =  0,  which  corresponds  to  diffusion  in  the 
direction  of  one  of  the  coordinate  axes.  To  show  that  (52)  is  only  conditionally  stable  for  certain  choices 
of  0  we  consider  the  special  case  0  =  ”,  so  we  require 

(«1  +«2)4  <  2 

4  +  a4  +  a4  +  \a\a\  ~ 

This  is  not  satisfied,  for  example,  when  oq  =  a2  =  1.  The  instability  is  easily  demonstrated  in  numerical 
simulations  and  avoiding  it  requires  using  time  steps  on  the  same  order  of  magnitude  as  an  explicit  scheme, 
k  =  0(h4). 
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